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PADE APPROXIMANT CALCULATION 
OF THE SINGULARITY IN THE 
MAGNETIC SUSCEPTIBILITY 
OF AN ISING SQUARE LATTICE 


by 

A. V. M. Ferris- Prabhu 
Goddard Space Flight Center 

SUMMARY 

The Pade approximant method was applied to Clapp's spin- 
spin correlation formulation for the magnetic susceptibility in an 
Ising square lattice, in order to investigate the nature of the 
critical point. For the exact cluster formulation it is found that 
as the cluster size increases, the critical point, namely the point 
at which the susceptibility has a singularity, approaches Onsager’s 
exact value of 0.414, and that x(t) ~ l/( t ~ t c ) - For the self- 
consistent formulation, twenty Pade approximants were calcu- 
lated. The value of the critical point converged very rapidly to 
0.41862770 and it was found that x(t) ~ l/(t - t c ) . 
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PADE APPROXIMANT CALCULATION 
OF THE SINGULARITY IN THE 
MAGNETIC SUSCEPTIBILITY 
OF AN ISING SQUARE LATTICE* 


by 

A. V. M. Ferris- Prabhu 
Goddard Space Flight Center 


INTRODUCTION 

In this work the Pade approximant method is applied to the formula given by Clapp* for the 
magnetic susceptibility of an Ising square lattice, in order to determine the location and nature of 
the critical point. Utilizing spin-spin correlation functions and on the basis of a cluster theory, 
Clapp has developed an exact method and a self-consistent method for performing the cluster cal- 
culations. For a detailed exposition on the derivation of these methods and their applications in de- 
scribing the behavior of magnetic systems, reference may be made to his thesis (Reference 1). 


THE MAGNETIC SUSCEPTIBILITY 

Fisher (Reference 2) gives the following high temperature formula for the susceptibility: 





n = 0 


T > T 


( 1 ) 


where N is the number of spins in the system, y the magnetic moment of each spin, k the Boltz- 
mann constant, T the absolute temperature, o- = 2S where S = 1/2 for the systems considered 
here, and (o- 0 z cr z y is the correlation function for the z components of the spins. Equation 1 is 
rewritten in the form (Reference 1): 


kT _ 1 3 + 4t 2 - 4t 4 

IV 2 ^ * 1 - : 


HO < l . 


( 2 ) 


♦This work was performed during a post doctoral appointment at Massachusetts Institute of Technology; and a preliminary report 
will appear in Physics Letters, 15(2): March 15, 1965. 

tSee also P. C. Clapp, "Self-Consistent Correlation Theory of ferromagnetism," Phys. Rev, 137(4A):A1295 _ Al305, 1965, which appeared 
while this work was in the process of being published. 
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where 


X(t) = 2<01>/(l + 2t 2 - 2t 4 ) , 

( 3 ) 

<01> = (01)/Z , 

( 4 ) 


t = tan h(/3j/2) , jB = l/kT, and J is the exchange integral. 

The correlation function (01) is a ratio of two polynomials in t, (01) and the partition function 
Z, whose degree depends on the number of spins in the cluster. (01) and z have been tabulated by 
Clapp for a number of two dimensional clusters. 


PADE APPROXIMANTS 

In solving physical problems, the answer often can be obtained completely in principle, but 
sometimes it cannot be expressed in a closed form. For instance, in the three dimensional Ising 
model, any finite number of terms in the power series expansion of the partition function can be 
computed. Knowing every term of the series is equivalent to knowing the function everywhere, 
and by analytic continuation it can be extended outside the region of convergence of the power 
series representation. However, in some cases it is unnecessary to compute the value of the 
function and then to continue it analytically. 

What is done is to replace the infinite series (of which only the first few terms may be known) 
by a ratio of polynomials. Such a ratio is called a Pade approximant (Reference 3) and the zeros 
of the numerator and denominator represent the zeros and poles of the function which is being 
approximated. Clearly, the more terms known in the series representation, the higher will be the 
degree of the Pade approximants that can be computed, and the more closely will the "true" func- 
tion be represented. However, the question "When is the Pade approximant method applicable?" 
is difficult to answer (Reference 4). 

An (N, M) Pade approximant, P (M) (x)/Q (N> (x) > has polynomials of degree M and N in the numer- 
ator and denominator respectively. In general it is sufficient to deal with the (N, N) approximants 
as they usually converge very rapidly (Reference 3). It will be seen later, that the degree of the 
numerator and of the denominator together cannot exceed the degree of the last term known in the 
series from which the approximants are being computed. 

The Pade method enables the real pole of a function to be found even if it lies outside the cir- 
cle of convergence of the power series from which the approximant was obtained. For example, if 
a function which has a known power series representation has poles in the complex plane inside a 
certain region, then by taking enough terms to represent these poles on the boundary of the circle 
of convergence, the Pade approximants can be computed. The real and complex poles of the func- 
tion outside this region are given by the zeros of the Pade denominator. Baker (Reference 3) ex- 
plains this in detail and reference may be made to his work. 
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A further advantage of the Pade method is that it enables the power with which a function di- 
verges to be found. This is evident from the nature of the Pade polynomials. But what if, as in 
the case of the susceptibility, which near the critical point behaves as i/(t - t c ) a , the pole is 
not a simple pole? In this case, forming the logarithmic derivative gives a function which has a 
simple pole: 




d _ -a 

at in * - (t-t c ) 


The power with which the susceptibility diverges at the critical point is obtained by evaluating 
the residue of the logarithmic derivative at the pole. For example, given the first few terms in 
the series expansion of the logarithmic derivative of the susceptibility, 




a 0 + a x t + a 2 t 2 + a 3 t 3 + * • • + a M t M + * • • , 


let us try to locate the critical point t c and a, where it is assumed that 


x(t) ~ i/(t - t c ) a . 

The series is first replaced by Pade approximants of the form 


A< N > (t) 
B< n > (t) 


N 



n "0 


N 



n= 0 


N = 1, 2, 3, • , M/2 . 


The critical point in the Nth Pade approximation is the least real positive value of t for which 
B (N) (t) becomes zero. The value of a can then be found by calculating the residue of the approx- 
imant at the critical point. 


-a 




AW (t) 
B< n > (t) 


This will become clearer by taking an example from Baker (Reference 3), who gives 


0(t) 


dt ln x(t) 


3 + 3t + 9t 2 + 15t 3 + 33t 4 + ••• , 
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for the logarithmic derivative of the susceptibility of a plane hexagonal lattice. Let 
B< 1 > (t) = b/ 1 ) + b/ x > t , 

A< l >(t) = a/ 1 ) + a/ x > t , 

®(t) B (1) (t) = 3b/ 1 ) + (3b/ 1 ) +3b 1 < 1 >) t + (9b/ 1 ) + 3b/ 1 )) t 2 + ••• , 

And since B (1) (t) and A< x ) (t) are each of degree 1, the following equalities may be made (Ref- 
erence 4): 


a 0 (1) 

= 3b 0 (1) . 

(5) 

a/ 1 ) 

II 

cr 

+ 

0* 

(6) 

0 

= 9b/ 1 ) + 3b/ 1 ) . 

(7) 


Inserting b/ n in terms of b/ 1 ) from Equation 7 into Equation 6 and eliminating b/ X) , we ob- 
tain the first Pade approximant: 

A< x ) (t) _ 2 (t ~ 2 ) 

B' x )(t) ” t _I 


And the residue is 


\ 3 / B< x > (t) 3 

t=l/3 


It is seen that the first Pade approximant locates the Curie point at t c = 1/3 and indicates 
that the power with which the susceptibility diverges is 1/3. By calculating higher Pade approx- 
imants, these values are refined. For the example we have just described, the critical point oc- 
curs at t c = 0.577 and the susceptibility behaves as l/(t - t c ) 175 . 

EXACT CLUSTER CALCULATIONS 

By defining a reduced susceptibility, 


Xr ( t ) 


kT 

N fJ . 2 


xt~2 X(t) , 
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Equation 2 can be rewritten as 


where 


XrCO = F(t)/G(t) , (8) 

F(t) = A(t) • Z , (9) 

A(t) = (3 + 4t 2 - 4t 4 ) , (10) 

G(t) = C(t) • Z - 2(01) , (11) 


C(t) = (l + 2t 2 - 2t 4 ) . 


( 12 ) 


The IBM 7 09 was used to evaluate the poly- 
nomials F(t) and G(t) and to locate the zeros of 
G(t). The smallest real positive value of t for 
which G(t) is zero is identified as the critical 
point t c . Table 1 gives t c for different clusters 
of increasing size. 


Table 1 

Location of Singularity. 


Number of Spins 
in Cluster 

6 

9 

20 

25 

30 

tc 

0.521 

0.473 

0.442 

* 

0.431 


It is evident from these results that the *This value was obtained and later misplaced, 
susceptibility has a singularity at a point which 

is identified as the critical point, and that as the cluster size increases, the critical point occurs 
at lower values of t . It is possible that in a cluster of an infinite number of spins, the critical 
point would converge to Onsager T s exact value (Reference 5) of 0.414. 


The Pade method (Reference 3) is now applied to the problem of determining the power with 
which the susceptibility diverges. Now, if 


*(t> ~ i/( t -t c )*, 


then 


a - -(t-t c ) • d(ln x}/dt l t=t< _ . 

Hence, if the expression x(t) = F(t)/G(t) behaves as l/(t - t c ) a in the neighborhood of the 
critical point, the power a should be obtainable by evaluating the residue of the logarithmic deriva- 
tive of the susceptibility at the critical point. This was done on the IBM 709 and the residue was 
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found to be exactly -1 for each cluster. This means that whenever the susceptibility is calculated 
from a finite cluster in the exact manner of Clapp, it will appear as a ratio of polynomials and 
thus will behave as l/(t - t c ) . This can easily be proven, for let 


x(t) 


E(t) 

(t-0 * 


where 


E(t) 


F(t) 

G(t) 


and 


G(t) 


G(t) 

(t-O ' 


Then 


d( In x) _ r (t) _ 1 

dt - E(t)-( t _g 


and the residue is 


(t - t c ) • d( In x)/dt 


(*“*=) 


E'-Ct) 1 

. E(,) (<-<.) 


“I • 


SELF-CONSISTENT METHOD CALCULATIONS 

The self-consistent method (Reference 1) may be explained briefly by considering an example 
of a two dimensional cluster of six spins. The bonds in Figure 1 represent the exchange coupling 

between spins located on the numbered sites. Unlike 
the usual cluster theory where all exchange couplings 
have the constant value J, in the self-consistent 
method only the coupling between spins 1 and 4 is 
assumed to have the value J. All other couplings 
are assumed to have a temperature dependent value 
Ji (/?) which is to be determined by requiring that 
at all temperatures the correlations between spins 1 
and 4 and between spins 1 and 2 (or any outside 
nearest neighbor pair) be equal. 



Figure 1— Spin cluster used for 
two-dimensional lattices. 
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The physical reasoning for this is that in an infinite lattice of which this cluster is a part, a 
spin pair like 1 and 2 for instance, will have, in addition to the coupling paths in the cluster, other 
paths in the rest of the lattice. If this temperature dependent additional coupling be represented 
by J' (/3), then the effective coupling between this spin pair will be given by J x (J3), that is by the 
constant coupling J plus J' (/ 3 ). In the case of the pair 1-4, it is assumed that the coupling from 
the rest of the lattice is accounted for by the augmented coupling between spins on the perimeter 
of the cluster, thus enabling the coupling to be kept at the constant value J. 

The correlation function (01>, a function of J and j x , now appears as a ratio of polynomials 
in t and t 19 where The requirement that in an infinite crystal all nearest 

neighbor correlations must be equal, leads to the self-consistency relation (Reference 1): 


1 - Vl - 4t 2 1 

1 1 2t ; t < 2 * 


Thus in the self-consistent method calculation, the susceptibility 


x R (O 


1 3 + 4t 2 - 4t 4 
3 1 - k 


K < 1 , 


can be put into the form 


Xr (t) 


z>* 

ft 


(13) 


where 


= ft 


L 

n = 0 

r 

« = ft 


F t n 


n 


G t n 


(14) 


I 



= A(t)C(t) ^E n t" , 


^G n t" = C(t)^]E n t n " 2^D n t n 


A(t ) = (3 + 4t 2 - 4t 4 ) 

C(t) = (l + 2t 2 - 2t 4 ) 


I> tn = i + 2 t (I>" tn ) + (E 


IS tn 


1 - yi-4t 2 

2t 


The logarithmic derivative of the susceptibility is 


d(ln x) 
dt 


where 




n = 0 

Z]°- ,n 


Q„ = (n + l)Q n + 1 . n = 0,1,2, 
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The IBM 709 was used to generate 41 terms in the power series expansion of t t and to com- 
pute the first 40 coefficients D n , E n , F n , G n , Q n , Q n and T n . In Appendix A the series representation 
for the susceptibility and its logarithmic derivative are given. For comparison the poles of the 
susceptibility as given by Equation 14 were located. It was found that the only real value for which 
there is a pole is at t = 0.41862991. 

The first 20 Pade approximants were then computed and their poles located. The critical 
point is identified as the smallest positive real value of t for which the approximant has a pole. 

At this point t c , the residue was calculated and the results are tabulated in Table 2. 

The fifth approximant locates t c at 0.4151549 which is in excellent agreement with Onsager’s 
exact value of 0.414. The sixth approximant behaves strangely with t c at 3.4049675 and a residue 
of 33.771. This result seems spurious though it is not clear why it appears. The next fourteen 
approximants behave well and t c converges rapidly to 0.41862770 with a slight oscillation in the 
last two decimal places. The residue rapidly becomes -1. It is interesting to observe that direct 
calculation of the singularity in the susceptibility locates t c at 0.41862991. 

Table 2 

zv 

Results derived from d(lnx)/dt - ~ — 



N 

location of t 

C 

residue 

N 

location of t c 

residue 

1 

0.83333333 

- 4.61 

12 

0.41862779 

- 1.00000 

2 

0.8259 

- 3.59 

13 

0.41862779 

- 1.00000 

3 

0.34672120 

- 0.336211 

14 

0.41862778 

- 1.00000 

4 

0.4072426 

- 0.793024 


(0.21132568* 

-11.3767 x 10" 7 ) 

5 

0.4151549 

- 0.915533 

15 

0.41862864 

- 1.00000 

6 

3.4049675 

33.771 

16 

0.41862795 

- 1.00004 

7 

0.41884024 

- 1.00726 

17 

0.41862788 

- 1.00001 

8 

0.41859223 

- 0.998611 

18 

0.41862784 

- 1.00001 

9 

0.41862295 

- 0.999777 

19 

0.41862816 

- 1.00001 

10 

0.41862277 

- 0.999769 

20 

0.41862770 

- 1.00002 

11 

0.41862652 

- 0.999937 





♦The first real pole occurs at this value. 


DISCUSSION AND CONCLUSIONS 

For the exact cluster calculations the value of the critical point appeared to converge to the 
exact value as the size of the cluster increased, and the susceptibility behaves as l/(t - t c ) . 
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For the self-consistent calculations, twenty Pade approximants were calculated and the crit- 
ical point was found to converge rapidly to 0.41862770 which is in excellent agreement with the 
value 0.41862991 obtained by directly calculating the location of the singularity in the susceptibil- 
ity. Here, as well, the susceptibility behaves as l/(t - t c ) which is consistent with the result of 
the exact cluster calculation. 

Although the behavior of the susceptibility at the critical point is not predicted accurately by 
the self-consistent method of Clapp, the critical point itself is located with excellent accuracy. 
Insofar as this is obtained by using the smallest cluster, the self-consistent formulation seems a 
useful one. It is proposed to investigate, in a similar manner, the singularity in the susceptibility 
and in the specific heat of a Heisenberg lattice. 
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Appendix A 


SERIES REPRESENTATION FOR THE SELF-CONSISTENT C LUSTER METHOD 
SUSCEPTIBILITY OF AN ISI NG SQUARE LATTICE 

X(t) = 3 + 6t + 16t 2 + 32t 3 + 60t 4 + 156t s + 384t 6 + 936t 7 

+ 2208t 8 + 5224t 9 + 12496t 10 + 29984t u + 71677t 12 

+ 170960t 13 + 408064t 14 + 975136t ls + 2330368t 16 

+ 5566560t 17 + 13294784t 18 + 31756544t 19 

+ 7 58630 39 1 20 + 181224890t 21 + 432897020t 22 

+ 1034067100t 23 + 2470133200t 24 + 5900589300t 25 

+ 14095108000t 26 + 33.669717 x 10 9 t 27 + 80.428637t 28 

+ 19.212452 x 10 10 t 29 + 45.8939 12 x 10 10 t 30 + 10.96 2941 x 10 u t 31 

+ 26.187793 x 10 11 t 32 + 62.556262 x 10 u t 33 + 14.943172 x 10 12 t 34 

+ 35.695605 x 10 12 t 3S + 85.268115 x 10 12 t 36 + 20.368479 x 10 13 t 37 

+ 48.655348 x 10 13 t 38 + 11.622579 x 10 14 t 39 + 27.763514 x 10 14 t 40 

SERIES REPRESENTATION FOR THE LOGARITHMI C DERIVATIVE OF THE SUSCEPTIBILITY 

^lnx(t) = 2 + 6 . 666666t + 8t 2 + 70.111114t 3 + 92t 4 + 223.4074t s + 408t 6 

+ 863. 60496t 7 + 2600t 8 + 6447.7359t 9 + 14368t 10 
+ 33151. 827t n + 82448t 12 + 199871. 87t 13 + 471008. 05t 14 

+ 1115647. 7t is + 2682975. 8t 16 + 6431553t 17 + 15325437t 18 
+ 36547323t 19 + 87409346t 20 + 20894122t 21 

+ 49.885005 x 10 7 t 22 + 11.912655 x 10 8 t 23 + 28.464097 x 10 8 t 24 
+ 68.002609 x 10 8 t 25 + 16.242093 x 10 9 t 26 + 38.795986x 10 10 t 27 
+ 92.680261 x 10 10 t 28 + 22.139729 x 10 10 t 29 + 52.884728 x 10 10 t 30 
+ 12.632670 x 10 10 t 31 + 30.176868 x 10 11 t 32 + 72.085836 x 10 n t 33 
+ 17.219383 x 10 12 t 34 + 41.. 1327 29 x 10 12 t 3S + 98.256681 x 10 12 t 38 
+ 23.471224 x 10 13 t 37 + 56.066712 x 10 13 t 38 + 13.372948 x 10 14 t 39 
NASA -Langley, 1965 G-633 11 
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